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ABSTRACT 

An optimal algorithm for the in-flight calibration of spacecraft gyroscope systems is presented. 
Special consideration is given to the selection of the loss function weight matrix in situations in 
which the spacecraft attitude sensors provide significantly more accurate information in pitch and 
yaw than in roll, such as will be the case in the Hubble Space Telescope mission. The results of 
numerical tests that verify the accuracy of the algorithm are discussed. 


I. INTRODUCTION 

A spacecraft gyroscope system, sometimes refeiTed to as the rate gyro assembly (RGA), is used 
to measure angular rotation rates of the spacecraft. This is required for proper control of the 
spacecraft, particularly for the proper positioning of spacecraft sensors with respect to desired 
targets. We present here an algorithm for RGA calibration that was used for the High Energy 
Astronomy Observatory (HEAO) missions and enhanced for use in the Hubble Space Telescope 
(HST) mission. Section II of this article presents the basic algorithm; Section III, the statistical 
weighting scheme; Section IV, the results of numerical tests of the algorithm; and Section V, our 
conclusions. 
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II. BASIC RGA CALIBRATION ALGORITHM 


Any RGA must be composed of at least three gyroscopes whose axis directions taken together 
completely span the space of possible rotations (i.e., pitch, roll, and yaw). An RGA consisting of 

K gyros produces as output a response "vector" R = (r t , r 2 , ... , r K ) T , where r k is the response of 
the kth gyro. The response vector is translated into a measured angular velocity, Q M , of the space- 
craft (in the spacecraft frame) via the relation 

= G q R - Dq , (1) 

where G 0 is the RGA 3-by-K scale factor / alignment matrix, and D 0 is the RGA drift rate bias. 

If G 0 and D 0 deviate from their true values, because of either poor initial calibration or temporal 
changes of the RGA, then will deviate from the true angular rate, Q. The goal of the algorithm 

is to determine correction matrices M and d that may be applied to G 0 and D 0 so that a modified 


equation (1) will yield the true angular rate: 

G = MG 0 , (2a) 

D = M D 0 + d , (2b) 

Q = GR-D=A/£2 M -d. (2c) 

The angular rate deviation, co, between the measured and true rates is given by 

co = n M - Q = -mCl M + d, (3) 


where m = M - / , and I is the identity matrix. The algorithm will solve for m and d. The cor- 
rection matrices m and d are dimensioned 3-by-3 and 3, respectively; we emphasize this because it 
ties directly to the fact that the algorithm being developed here provides correction information for 
the RGA as a whole in its capacity as a device for measuring three-dimensional angular motion. 
Unless the RGA under consideration consists of only three gyros, m and d will not contain suffi- 
cient information to allow separate calibration updates of scale, alignment, and drift for the indi- 
vidual gyros. 
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We wish now to relate the angular rate deviation as integrated over some calibration maneuver to 
the difference in attitude changes as determined for the maneuver by the RGA and by some inde- 
pendent attitude sensing devices (e.g., fixed-head star trackers). This derivation is conveniently 
done using quaternion notation; see Reference 1 for a discussion of the mathematics of quaternions. 
First, an expression for the time derivative of a quaternion is required. If Q(t) is the quaternion 
representing spacecraft attitude at time t, then the quaternion representing a change in attitude over a 
time interval At is given by Q -1 (t) Q(t+At). If At is small, the attitude change may be expressed as 

Q'Ht) Q(t+At) = 1 + q(fiAt/2), ** (4) 

where 1 is the identity quaternion, and q(QAt/2) is a differential quaternion with vector component 
QAt/2 and scalar component zero. Defining AQ(t) as Q(t+At) - Q(t) and combining this defini- 
tion with equation (4) yields 

AQ(t) = Q(t) q(QAt/2) . (5) 

Dividing equation (5) by At produces the desired quaternion time derivative, Q'(t): 

Q'(t) = AQ(t)/At = Q(t) q(Q/2) . (6) 

Equation (6) applies as well for the quaternion time derivative corresponding to the attitudes as 
measured by the RGA, with subscript M placed appropriately. 

Next, the time derivatives specified above are used to construct the time derivative of the attitude 
error quaternion and the definite integral of that quaternion over the time of the maneuver. The 
attitude error quaternion, 5Q, is defined as 

SQ = Q m < <V Q > <V = Q Q m ‘ • O) 

which is a quaternion expressing a rotation from the RGA-determined postmaneuver attitude to the 
true postmaneuver attitude, transformed to the premanuever reference frame. It follows by the 
chain rule of differentiation that 

5Q’ = QQm' 1 + QQ’m 1 <&0 

- Q q(fi/2) Q M 1 + Q q '<n M /2) Q m -i . (8b) 
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Combining equation (8b) with the relations CO = Q M -Q and q _1 (0 M /2) = q(-^l M /2) produces 

5Q’ = Q q(-co/2) Q M i . (9) 

Integrating both sides of equation (9) over the maneuver yields 

6Q-1 = J Q q(-co/2) Q M _1 dt , (10) 

where the constant of integration (i.e., the identity quaternion) is removed from 8Q because the 
integral in equation (10) is a definite integral. The attitude error quaternion can be expressed in 
terms of the rotation Q R1 _1 Qj^ from the first reference attitude to the second (i.e., as determined 

using the attitude sensors against which the RGA is being calibrated) and the rotation Qqj ^02 
between the first RGA-propagated attitude and the second. Equation (10) thereby becomes 

(Qri- 1 Qr 2 ) (QG2- 1 Qgi> - 1 = Jqq(-co/2)Q m -i dt . (11) 

No approximations have been made in the derivation to this point. We now make two approxi- 
mations, each of which is accurate to first order in the error. First, Q M is substituted for Q in the 
integrand in equation (11). This substitution yields 

(Qri" 1 Qr2^ (Qg2* Qgi) " 1 = J Qm Q( _ gV2) Qm’ 1 dt • (12) 

The integrand in equation (12) is simply the quaternion representation for a rotation of the vector 
-co/2 through a rotation defined by Q M , i.e., the rotation that transforms -co/2 from spacecraft 

coordinates at time t to spacecraft coordinates in the premaneuver reference frame. Equation (12) 
can therefore be written in matrix notation as 

Zj = - 1/2 J 7 1 . co dt (13a) 

= 1/2 \t. (m Q m - d)dt, (13b) 

where Z ; is the vector component of 8Q, T- x is the matrix for transforming vectors to premaneuver 
spacecraft coordinates, and i is a subscript designating maneuver number. The second 
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approximation is made implicitly in the definition of as the vector component of 8Q { ; the fourth 
component of 5Qj, which is actually equal to the cosine of the error rotation angle, is approximated 
as equal to 1. Because of the two approximations made in going from equation (1 1) to (13b), the 
calibration algorithm described here will be inherendy iterative; the vector Z ; and matrix T { must be 
reevaluated on each iteration. 


Equation (13b) is linear in the unknowns m and d and thus lends itself naturally to standard 
least- squares techniques. First, the matrix equation that represents equation (13b) applied to N 
calibration maneuvers is written as 

Z = H x , (14) 

where Z and the state vector x are defined via 

Z = {Z 1 t,Z 2 T,...,Z n T} t , (15) 

x — 1/2 { m 11? m 12 , m^, m^, ^23* ^ 31 * ^ 32 * ^33* ^ 2 * ^3 ^ * (16) 

and H is a 3N-by-12 matrix of the form 


U 1 

H ~ u N -r N ' 

Each U i is a 3-by-9 matrix, the components of which are given by 
(Uj k+3(l-l) = ^ k^i dt » 


(17) 


(18) 
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(19) 


and each Y- t is a 3-by-3 matrix with components 

< Y jk>i = /(T jk )i *. 

The least-squares solution for the state vector x in equation (14) minimizes the linear Bayesian 
weighted least-squares loss function, J, given by 

J=1/2[E t WE + ( x - x a ) T S a ( x - x a ) ] , (20) 

where 

E = Z-Hx, (21) 

W and 5 a are symmetric nonnegative definite weighting matrices, and x a is an a priori estimate 
of x. The desired solution for x is given by 

x = (^WH + SJ- 1 (H t WZ + S a x a ). (22) 

Equations (20) and (22) include S a and x a for mathematical completeness. In what follows, we 

will assume S a = [0], i.e., no a priori knowledge of x. In applying equation (22) to determine x, 

it is clear that at least four calibration intervals are required and at least three of these must span the 
space of possible rotations. An acceptable minimum set of calibration intervals would be one 
maneuver each of pitch, roll, and yaw, together with a period of constant attitude to define the drift 
rate bias. In selecting calibration maneuvers to be used for the algorithm, a user should be aware 
that, at least in the equations specified above, a rotation of greater than 180 degrees is indistin- 
guishable from a smaller rotation in the opposite direction. The use of such large rotations could 
lead to errors in the calibration and should therefore be avoided. The basic algorithm being dis- 
cussed here lends itself easily to being broken into two separate algorithms, one to determine the 
scale factor / alignment portion of x and a separate one to determine the drift rate bias. 
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III. SPECIFICATION OF THE WEIGHT MATRIX 


In principle, the specification of the weight matrix W in equation (22) depends on the scale size 
of random errors associated with the RGA itself, as well as the errors associated with the determi- 
nation of reference attitudes. In practice, random errors associated with the RGA tend to be much 
smaller than those of the reference attitude sensors. This is true in particular for the sensors used 
for the HST mission. The HST uses two types of sensors for high-accuracy attitude determination: 
fixed-head star trackers (FHSTs) and fine guidance sensors (FGSs). The three FGS fields of view 
are clustered tightly (within 14 arc-minutes) about the principle axis of the spacecraft (hereafter 
called the Vl-axis). When calibrated, the accuracy of the FGSs should be better than 0.010 arc- 
second. In most circumstances, however, their effective accuracy will be limited by the accuracy of 
the reference star catalog against which the FGS observations are compared during attitude compu- 
tation; this will be about 0.3 arc-second. Because of the tight clustering of the FGSs about the 
Vl-axis, this accuracy of 0.3 arc-second pertains only to the pitch and yaw components. The roll 
accuracy is determined by the accuracy of the FHSTs, whose fields of view are more widely 
distributed in direction about the spacecraft than those of the FGSs. The accuracy of the FHSTs 
used for the HST is about 10 arc-seconds. The fact that the reference attitudes for the HST are 
substantially better in pitch and yaw than in roll presented a special problem for the HST ground 
software. For spacecraft that have equal attitude accuracies about all axes, setting the weight matrix 
in equations (20) and (22) to the identity matrix, i.e., treating the accuracies of all components of all 
maneuvers as equal and uncorrelated, would be legitimate. This is not the case for the HST; conse- 
quently, significantly improved results for HST RGA calibration can be expected if a proper weight 
matrix is used. The need for a proper weight matrix is enhanced by the possible requirement to 
combine data sets for maneuvers in which some of the attitudes were determined using only FHSTs 
or only FGSs. This possibility arises because of potential sensor occultation by the Earth during 
parts of the spacecraft orbit 

The net effect of the considerations discussed above is that the simple product E T E (E defined in 
equation (21) ) does not represent a squared sum of normalized, independent random variables as is 
required for an optimized least-squares loss function (e.g., see Reference 2). The determination of 
W depends upon the measurement uncertainties of the components of E, both in magnitude for the 
individual components and in any correlation of errors between the individual components. The 
3N- vector E is composed of N 3-vectors, the relation being 
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(23) 


E = (ej T , e 2 T , ... , e N T ) T , 

where e ; is the error associated with the ith calibration maneuver. If die errors for the maneuvers 
are assumed to be independent, then W can be represented as a block diagonal matrix, with each 
block being a 3-by-3 matrix, w^. The loss function becomes 


J = 2^, (24a) 

jj = 1/2 eT h>j e ; . (24b) 

The assumption is not strictly valid in cases where the same attitude data are used at both the end 
of one maneuver and the beginning of another. In such cases W would appropriately contain 
elements representing a cross-correlation of errors between maneuvers. Because this complication 
is both awkward to include computationally and of substantially smaller effect than that of the 
accuracy asymmetries in attitude produced by the use of FHSTs and FGSs together, we will 
neglect it. We will also make use of the fact that the random RGA errors are negligible when 
compared with the reference attitude errors. For notation simplification, hereafter we will suppress 
the subscript i (for maneuver number) unless it is explicitly required. 

Using the approximations described above, the components of e may be written as 


e j = a j + Et jk b k> ( 25 > 

where a^ is the premaneuver attitude determination error about the jth spacecraft axis, b k is the 
postmaneuver error about the kth axis, and is t^ the (j,k) component of the postmaneuver-attitude- 
to-premaneuver-attitude transformation matrix. The numbers a p a 2 , a 3 , bj, b 2 , and b 3 may be 

assumed to be independent random variables, but the numbers e 1? e 2 , and e 3 will in general be cor- 
related because of the mixing of the postmaneuver errors via the maneuver transformation matrix. 
For cases where the and b k are all approximately equal (as would be true for attitudes determined 
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using a number of well-separated sensors of equal accuracy), the correlation between the compo- 
nents of e would be small because of the combining in each of a number of random variables in 

different ways. Here, however, we are particularly interested in contexts where a t and b } (the roll 

errors) are large compared with a^ a 3 , b 2 , and b 3 . For such cases a maneuver with significant 
components about all axes will have errors given to first order by 


til b l 

+ 

0(8) , 

(26a) 

hi b i 

+ 

0(8) , 

(26b) 

hi b i 

+ 

0(8) , 

(26c) 


where 8 is a generic random variable with variance like that of a 2 , b 2 , or b 3 . The correlation, 

particularly between e 2 and e 3 , is clear. The goal now is to construct from the components of e 

three independent, normalized random variables that may be used in defining a least-squares loss 
function. The components of the weight matrix w can then be solved for by setting this new loss 
function equal to 1/2 (e T w e). For this derivation we will consider two cases: (1) where the 

maneuver has a significant nonroll component (i.e., t 2 i 2 a bl 2 > o g 2 and/or t 31 2 cr bl 2 > cr 5 2 , 
where G bl 2 and c g 2 represent the variances of bj and 8, respectively) and (2) where the maneuver 
is essentially pure roll (i.e., t^c^ 2 < o g 2 and t 31 2 G bl 2 < a g 2 ). 


In case 1 (a maneuver with a significant nonroll component), the quantities e^ e 2 and e 3 may be 

used to construct three essentially independent random variables <a t >, <bj>, and c, where 
<a t > and <bj> are estimators of a 1 and bj, and c is a variable with variance like that of 8. We 
specifically construct c to eliminate the large roll errors: 

C = *31 e 2 ' *21 e 3 * ( 27 ) 

The variance of c over an ensemble of maneuvers is given by 

o c 2 = var { t 31 e 2 - e 3 } (28a) 

= (t 31 2 C 2 + t 21 2 C 3 )o g 2 , (28b) 
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where 

C 2 = 1 + hi + t 23 2 ’ ( 29a ) 

C 3 = 1 + t 32 2 + t 33 2 . (29b) 

We construct <bj> to be a good estimator of b t : 

<bj> = [ t^ 2 / (t 2 i 2 + ^j 2 ) ] (e 2 / 1-21 ) + [ t^ 2 / (^i 2 + tjj 2 ) ] (e 3 / tjj ) . (30) 

The quantities e 2 /t 2l and e 3 /t 31 estimate bj to within S/t^ and 5/tjj, respectively. The 
quantities in brackets are normalized weights inversely proportional to the square of the uncertainty 
of the corresponding estimate. For convenience we rewrite <b t > as 


<bj> = B 2 e 2 + B 3 e 3 , 


(31) 


where 

B 2 = *21 / ( * 2 i 2 + hi 2 ) 5 (32a) 

B 3 = *31 / ( hi 2 + hi 2 ) • (32b) 


The variance of <bj> over an ensemble of maneuvers can be shown to be 
°<W> 2 = V + (V<1 + B 3 2C 3)° 6 2 - 


(33) 


The variance of <b t > increases as the maneuver approaches being pure roll; if t 21 2 °bi 2 = 
t 3 i 2 a b i 2 = cj 6 2 , then a <bl> 2 = ( 1 + C 2 /4 + C 3 /4 ) o M 2 = 2a bl 2 . Finally, we construct 

<a : > to be a good estimator of a t : 


<a t > — ej * tjj <b t > — ej — A 2 e 2 — A 3 ^ 3 j 


where 


A? - *11 B 2 ’ 


A 3 - *n B 3- 


(34) 


(35a) 

(35b) 
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( 36 ) 


The variance of <aj> over an ensemble of maneuvers can be shown to be 

°<al> 2 = V + <*12 2 + l 13 2 + A 22 ^ + A 32 C3 ) . 

The variance of <aj> also increases as the maneuver approaches being pure roll; if t 21 2 o bl 2 = 
hi^bi 2 = a 5 2 ’ then a <ai > 2 = a ai 2 + ( C 2^ 4 + C 3 /4 ) o M 2 = a ^ + a bl 2 Although 
< aj > is defined using <bj>, it is specifically tailored to remove the correlation with b t from e r 
To lowest order <a 1 > = a 1 , <bj> = bj, and c = f^, a 3 , b 2 , b 3 ), from which it is clear that 
<aj>, <bj>, and c are essentially independent. To find expressions for the weight matrices w, 

we construct a loss function from the squared sum of <a,>, <bj>, and c after normalization and 
set it equal to the original loss function, i.e.. 


2 j = e T w» e = c/ct c 2 + <b 1 > 2 /a <M> 2 + <a x > 2 / or <al> 2 . (37) 

The corresponding elements of w are 

W 11 = 1 / a <al> 2 ’ ( 38a ) 

W 22 = A 2 2 / a <al> 2 +B 2 2 /°<bl> 2 + l 31 2 / a c 2 ’ ( 38b ) 

w 33 = A 3 2 /a <al> 2 +B 3 2 /a <bI> 2 + t^ 2 / a 2 , (38c) 

w i2 = w 2i = ( 38d ) 

w 13 = W 31 = - A 3 2 / a <al> 2 > ( 38e ) 

W 23 = W 32 = A 2 A 3/ a <al> 2 + B 2V a <bl> 2 + *21 *31 ' G c • ( 38f ) 

In case 2 (an essentially pure roll maneuver), the components of e can be expressed as 

ej = a^ + bj , (39a) 

e 2 = a 2 + t 21 b 1 + (cos 0 ) b 2 + (sin 0 ) b 3 , (39b) 

e 3 = a 3 + t 31 bj - (sin 0)b 2 + (cos0)b 3 , (39c) 
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where 0 is the roll angle. Because b 2 and b 3 are assumed to have equal variance, the variables 
r + and r_, defined as 

r + = (cos 0 ) b 2 + (sin 0 ) b 3 , 
r. = - (sin 0 ) b 2 + (cos 0 ) b 3 , 

are independent random variables with the same variance a g 2 . 
of the form 


e l = a l 

+ 

b l’ 

(41a) 

e 2 = r 2 

+ 

*21 b l ’ 

(41b) 

e 3 = r 3 

+ 

t 3] bj , 

(41c) 


where r 2 and r 3 are independent random variables with variance 2o g 2 . The components of e 

are mildly correlated via bj; bj contributes half of the variance of ej and at most one- third of the 
variances of e 2 and e 3 (for t^ 2 = t^ 2 = G s 2 / G bl 2 ). We neglect this mild correlation for 
maneuvers that are essentially pure roll by treating them as exactly pure roll, i.e., by setting 
t 2 i = t 3 j = 0. The weight matrix elements that follow from this assumption are 


w n = 

l/( 

<j 2 -f (y 2 ) 

(42a) 

W 22 = 

W 33 

= 1/2g 5 2 , 

(42b) 

W 12 = 

W 21 

= w 23 = w 32 = w 13 = w 31 = 0. 

(42c) 


The equations specified above provide the functional relationship between the elements of w and 
the uncertainties in attitude determination with respect to the spacecraft axes. These uncertainties 
can be derived from the attitude covariance matrix (e.g., see Reference 3), given by 

P = G 2 [ / - £ ( G 2 / G k 2 ) V k VJ ] 4 (43) 

k=l,n 


(40a) 

(40b) 

The components of e are therefore 
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where 


n = number of measurements, 


V k = kth star vector used for attitude determination, as expressed in the spacecraft frame, 
a k = uncertainty of kth measurement, 


CT t 2 = [ 2 ( C7 k 2 r 1 ]- 1 . 


For spacecraft such as the HST, whose sensor orientations allow attitude determinations with 
uncorrelated estimates of pitch, roll, and yaw, the diagonal elements of the matrix P may be used as 
the attitude error variance required in the expressions for the elements of w. 


IV. NUMERICAL TESTING 

The scheme specified above for applying statistical weights to RGA maneuver data has been 
implemented in the Payload Operations Control Center (POCC) Applications Software Support 
(PASS) system to be used in support of the HST mission. As indicated at the end of Section II, the 
algorithm was implemented in a way that allows independent calibration of the RGA scale factor / 
alignment and the RGA drift rate bias. RGA, FGS, and FHST data appropriate for a number of 
different calibration intervals were generated using the PASS attitude simulator (this is a simulator 
that produces data like that expected from the HST). The FGS and FHST errors were of order 0.5 
and 10 arc- seconds, respectively. The data consisted of nine independent 90-degree maneuvers 
(three each of pitch, roll, and yaw) and one extended period of constant attitude data. The 
90-degree maneuvers would require approximately 16 minutes of spacecraft time, whereas the 
constant attitude data represented approximately 45 minutes of spacecraft time. These data were in 
turn processed through the PASS attitude determination and RGA calibration software. Final 
residuals were calculated for each maneuver component by comparing the RGA-measured 
maneuver with the "true" maneuver as supplied to the simulator. The results were good. The final 
residuals for the components of E were in magnitude appropriate for the sensor type governing 

those residuals (e.g., e 2 for a pure pitch maneuver was in size like the simulated FGS errors, 
whereas c 1 was like the simulated FHST errors). Furthermore, the scheme described above allows 
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for the combining of data sets in which different attitudes were determined with various different 
sensor combinations (i.e., FHST and FGS, FHST only, and FGS only). To test this aspect, we 
processed our simulated data through two RGA calibration scenarios: (1) with both FHST and 
FGS data used for all attitude determinations and (2) with each maneuver processed twice, once 
with FHSTs only and once with FGSs only. The final accuracy of the RGA calibration was 
essentially the same for both scenarios. For the nine-maneuver simulation, the largest deviation 
(when comparing the measured maneuver magnitude after calibration with the true maneuver 
magnitude) was found to be about 30 parts per million. This is both appropriate for the magnitude 
of the sensor errors and sufficiently accurate to support the needs of the HST mission. 


V. CONCLUSIONS 

We have presented a general algorithm for the calibration of a spacecraft rate gyro assembly, as 
well as a data weighting scheme that produces a statistically optimal solution. The weighting 
scheme, although explicitly tailored for use during the Hubble Space Telescope mission, is 
applicable to any three-axis stabilized spacecraft. Numerical simulations demonstrate that the 
algorithm works as expected in theory and is capable of supporting the needs of the HST mission. 

The work reported in this article was supported in part by NASA contract NAS 5-26685 for the 
development of ground support software for the Hubble Space Telescope mission. 
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